# 假定属性个数K=3 G=5 首先将作答数据模拟 基于DINA模型
rm(list = ls())
library("psych")
library("poLCA")
library("ggplot2")
library("dplyr")
library("arrangements")
library("gtools")
set.seed(555)
N <- 100  #被试数
K <- 3              #属性个数
J <- 15             #题目数量
source(file = "functions.R") #加载自定义函数
#---------------------------------------------------------------------------------------

#题目参数（0.1-0.2），DINA模型
S <- runif(15, min = 0.1, max = 0.2) #失误参数
G <- runif(15, min = 0.1, max = 0.2) #猜测参数
ITEM_P <- cbind(S,G)
#真值题目的Q可能情况（应该是将属性关系图作为线索？应该借助贝叶斯网络结构来表示BN）
#Attribute hierarchy relationship 分别为属性1和属性3扩散
# # 邻接矩阵是否可以通过属性层级关系得出
#真实邻接矩阵得到可达矩阵和属性掌握模式 010 110
A <- matrix(c(0,1,1,
              0,0,0,
              0,0,0),nrow = 3,byrow = T)
# B <- matrix(c(0,0,1,
#               0,0,1,
#               0,0,0),nrow = 3,byrow = T)
l1 <- ad_Reachability(A) #得到可达矩阵（可自达）和属性掌握模式
# l2 <- attribute_def(B)
reach_ability_A  <- l1$reachability_matrix
attribute_A <- l1$attribute_matrix #属性掌握模式
item_A <- l1$test_item
# reach_ability_B  <- l2$expanded_matrix
# attribute_B <- l2$attribute_matrix #属性掌握模式
# 抽取测试题目 真实场景下上述信息都是未知

sampled_items <- lapply(J, function(size) {
  item_A[sample(1:nrow(item_A), size, replace = TRUE), ,drop = F ]
})    #根据设定的Q进行可放回随机抽取出Q矩阵

#Q[sample(1:nrow(Q), 15, replace = TRUE), ,drop = F ]

#把前3行替换为可达矩阵 保证可识别性
sampled_items[[1]][1:3,] <- t(reach_ability_A)
rownames(sampled_items[[1]]) <- paste0("Q",1:J)
# 从设定的属性掌握模式中抽取被试，共100人（多个提取）
sampled_people <- lapply(N, function(size) {
  attribute_A[sample(1:nrow(attribute_A), size, replace = TRUE), ,drop = F ]
})
# 单个sampled_people <- attribute[sample(1:nrow(attribute), N, replace = TRUE), ,drop = F ]
# 将行索引作为每次的限速进行抽取

# 人不变 题变即Q矩阵进行变化 方案一 前三行替换为可达 其他暴力替换101变010
# Q2 <- as.data.frame(sampled_items[[1]])
# Q2 <- t(apply(Q2, 1, function(row) {
#   if (all(row == c(1, 0, 1))) {
#     return(c(0, 1, 0))  # 替换为 c(0, 1, 0)
#   }
#   return(row)  # 否则返回原行
# }))
# # 再将可达矩阵放入
# Q2[1:3,] <- t(reach_ability_B)

#-----------------------------------------------------------------------------------------------
#作答矩阵的生成 共100个人在15道题目上基于DINA模型生成作答,先得出概率再进行判断
# pro <- DINA_model(sampled_people[[1]],sampled_items[[1]],S,G)
# # 生成作答
# Response <- matrix(NA, nrow = N, ncol = J)
# for (i in 1:N) {
#   for (j in 1:J) {
#     Response[i, j] <- ifelse(pro[i, j] > runif(1), yes = 1, no = 0)
#   }
# }
# 得到作答数据，先用平行分析进行估计属性个数（后期要进行优化！！）
Response <- DINA_model(attribute_profile = sampled_people[[1]],
                       Q_matrix = sampled_items[[1]],S,G)


fa_factors  <- fa.parallel(Response,show.legend = F,plot = F)$nfact

#在不存在任何属性关系的情况下，Q的情况 自定义函数
Q_est_num <- (2**fa_factors) - 1  #完备矩阵的个数
# Q_est <- generate_complete_Q_matrix(K)[-1,] #完备Q矩阵，没有层级关系
# people_est <- generate_complete_Q_matrix(K) #完备Q矩阵，没有层级关系

# LCA对人进行分类，得到属性掌握模式的类别数class，response_Initial为0，1作答的
response_Initial <- as.data.frame(Response)  #原始作答矩阵存储
Response[Response == 1] <- 2   #数据中不能出现0，进行数据转换将0变成1，1变成2
Response[Response == 0] <- 1
colnames(Response) <- paste0("Q",1:J)
# 处理成LCA公式输入格式
f1  <- as.formula(paste("cbind(", paste(colnames(Response), collapse = ", "), ") ~ 1"))

# 如果只单独考虑BIC的指标进行LCA的分析（最好考虑交叉验证最佳类别，对人能力进行分类）
# min_bic <- 100000  # 最小标准的信息量（越大越好，会被算法替代）
# for (i in 2:Q_est_num) {
#   lc <- poLCA(f1, data = as.data.frame(Response), nclass = i, maxiter = 1000,
#           tol = 1e-10, nrep = 10, verbose = FALSE, calc.se = TRUE,na.rm = FALSE)
#   if (lc$bic < min_bic) {
#     min_bic <- lc$bic
#     LCA_best_model <- lc
#   }
# }


# LCA找到最适合的类别
max_II <- -100000  # 达到收敛的公差值
min_bic <- 100000  # 最小标准的信息量
LCA_best_model <- NULL
LCA_resulst <- vector(mode = "list", length = length(2:Q_est_num))

#BIC作为主要的类别选择指标 
for (i in 2:Q_est_num) {
  lc <- tryCatch({
    poLCA(f1, as.data.frame(Response), nclass = i, maxiter = 3000,
          tol = 1e-5, nrep = 20, verbose = FALSE, calc.se = TRUE)
  }, error = function(e) {
    # 如果模型拟合失败，跳过当前循环
    return(NULL)
  })
  LCA_resulst[[i - 1]] <- lc   # 如果i为2那就是[[1]]，拟合3就是[[2]]
  if (!is.null(lc) && lc$bic < min_bic) {
    min_bic <- lc$bic
    LCA_best_model <- lc
  }
}
people_best_class <- length(LCA_best_model$P)#结果没有准确，考虑采用区间范围G
# 如果找到了最佳模型 准确的话应该是 5
LCA_result_people_mistake <- vector(mode = "list", length = length(1:3))
# people的分类 每一类的人作为真值与之后的做题的分类做比较
for (i in 1:3) {
  LCA_result_people_mistake[[i]] <- which(LCA_best_model$predclass == i)
}
# 错误识别 110 100 000的 可以把识别出来的提出来，再放进去识别？
#----------------------------------------------------LCA分类的结果---------------------------

# 需要的类别分类结果 即将人分为5类
need_LCA_result <- LCA_resulst[[4]][["predclass"]]
LCA_result_people <- vector(mode = "list", length = length(1:5))
# people的分类 每一类的人作为真值与之后的做题的分类做比较
for (i in 1:5) {
  LCA_result_people[[i]] <- which(need_LCA_result == i)
}
LCA_result_ture <- LCA_resulst[[4]]
# 提取每个类别的作答概率
category_probs <- LCA_result_ture$probs

# 获取每个类别的人群pr(1)答错，pr（2）是答对概率
#每个类别在各个题目上的答对概率
category_probs_class_ture <- matrix(data = 0,nrow = 15,ncol = 5)
for (i in 1:15) {
  category_probs_class_ture[i,] <- category_probs[[i]][,2]
}
category_probs_class_ture <- round(category_probs_class_ture,4)
row.names(category_probs_class_ture) <- paste0("Q",1:15)

# 画图
lcmodel <- reshape2::melt(category_probs, level = 2)#将 LCA 模型的概率矩阵重塑为长格式，指定合并的维度
lcmodel$L2 <- factor(lcmodel$L2, levels = paste("Q", 1:15, sep = ""))# 调整顺序
zp1 <- ggplot(lcmodel,aes(x = L2, y = value, fill = Var2)) #初始化 ggplot 对象，设置 x 轴为 L2（问卷项目），y 轴为 value（每个类别的概率），用 Var2 作为填充颜色。
zp1 <- zp1 + geom_bar(stat = "identity", position = "stack")# 加堆积条形图层，stat = "identity" 表示 y 轴的值直接使用数据中的值


zp1 <- zp1 + facet_grid(Var1 ~ .) #按 Var1 创建面板，便于比较不同类别的条形图
zp1 <- zp1 + scale_fill_brewer(type="seq", palette="Greys") +theme_bw()
zp1 <- zp1 + labs(x = "题目",y="作答概率情况", fill ="概率")
zp1 <- zp1 + theme( axis.text.y=element_blank(),
                    axis.ticks.y=element_blank(),
                    panel.grid.major.y=element_blank())#隐藏 y 轴的文本、刻度和主要网格线，简化图形
zp1 <- zp1 + guides(fill = guide_legend(reverse=TRUE))
print(zp1)

# 题目分类简易版本 只需要将作答放入得到分类 准确率挺高的
matched_groups <- match_rows_end(response_Initial,4)

# matched_groups <- match_rows_end_no_number(response_Initial)
matching_rows <- list()
# 答对各类题目的分类 例如答对第一类题目的人是哪些！
for (i in 1:length(matched_groups)) {
  people_lable <- matched_groups[[i]]
  response_Initial <- as.data.frame(response_Initial)
  colnames(response_Initial) <- paste0("Q",1:J)
  matching_rows[[i]] <- which(apply(response_Initial[, people_lable], 1, function(row) all(row == 1)))
}#提出分类，依据分类提出作答数据，对作答数据的每一行进行严格判定

#---------------------------------------------------------------分割线-----
HDD_factor <- HDD(response_Initial,IMP = attribute_A,Qmatrix = sampled_items[[1]])
result_1 <- distinguish_Q(LCA_result_people = LCA_result_people,HDD_factor)

# 穷举过程
linjie_matrix <- generate_linjie_matrix(fa_factors)
linjie_matrix <- generate_linjie_matrix_2(fa_factors)
att_matrix <- list()
ture_lj_list <- list()
index <- 0
for (i in 1:length(linjie_matrix)) {
  lj <- linjie_matrix[[i]]
  attribute_matrix <- ad_Reachability(lj)$attribute_matrix
  if (nrow(attribute_matrix) == length(LCA_result_people)) {
    index <- index + 1
    # 改
    att_matrix[[index]] <- attribute_matrix
    ture_lj_list[[index]] <- lj 
  }
}
#  错误LCA
for (i in 1:length(linjie_matrix)) {
  lj <- linjie_matrix[[i]]
  attribute_matrix <- ad_Reachability(lj)$attribute_matrix
  if (nrow(attribute_matrix) == length(LCA_result_people_mistake)) {
    index <- index + 1
    att_matrix[[index]] <- attribute_matrix
    ture_lj_list[[index]] <- lj 
  }
}


# 生成典型项目考核矩阵
zi_matrix_list <- lapply(ture_lj_list, function(lj) ad_Reachability(lj)$test_item)

# 计算可达矩阵
re_matrix <- lapply(ture_lj_list, function(lj) ad_Reachability(lj)$reachability_matrix)
# all_Q_matrix <- list()
# for (i in 1: length(zi_matrix_list)) {
#   one_Q <- genarate_Q_matrix(zi_matrix_list[[i]],matched_groups,re_matrix[[i]],
#                              number_items = J)
#   all_Q_matrix[[i]] <- one_Q
# }

all_Q_matrix <- list()
for (i in 1: length(zi_matrix_list)) {
  one_Q <- genarate_Q_matrix(zi_matrix_list[[i]],matched_groups,re_matrix[[i]],
                             number_items = J)
  all_Q_matrix <- append(all_Q_matrix,one_Q)
}

all_Q_lv <- list()
Q_lj_index <- c()   # 用于记录每个 Q 对应的 lj 编号
in_number <- 1

# for (i in 1:length(all_Q_matrix)) {
#   browser()
#   for (j in 1:length(all_Q_matrix[[i]])) {
#     Q_ma <- all_Q_matrix[[i]][j]
#     Q_ma <- do.call(rbind,Q_ma)
#     
#     HDD_factor <- HDD(response_Initial, att_matrix[[i]], Q_ma)
#     one_lv <- distinguish_Q(LCA_result_people = LCA_result_people, HDD_factor)
#     
#     all_Q_lv[[in_number]] <- one_lv
#     Q_lj_index[in_number] <- i  # 记录此结果属于第几个 lj
#     in_number <- in_number + 1
#   }
# }
for (i in 1:length(all_Q_matrix)) {
  browser()
  for (j in 5:6) {
    Q_ma <- all_Q_matrix[[i]][j]
    Q_ma <- do.call(rbind,Q_ma)
    
    HDD_factor <- HDD(response_Initial, att_matrix[[i]], Q_ma)
    one_lv <- distinguish_Q(LCA_result_people = LCA_result_people, HDD_factor)
    
    all_Q_lv[[in_number]] <- one_lv
    Q_lj_index[in_number] <- i  # 记录此结果属于第几个 lj
    in_number <- in_number + 1
  }
}



# 计算每组的最小重合度值
min_zhi <- sapply(all_Q_lv, min)

# 找出最高最小值对应的结果
max_weizhi <- which(min_zhi == max(min_zhi))
best_lj_index <- Q_lj_index[max_weizhi]
linjie_end <- ture_lj_list[unique(best_lj_index)]



















#真值
attribute_people_classify <- list()
for (i in 1:nrow(attribute_A)) {
  attribute_people_classify[[i]] <- which(apply(sampled_people[[1]], 1 ,function(row) all(row == attribute_A[i,])))
}
name_list <- apply(attribute_A,1,function(x) paste(x, collapse = "_"))
names(attribute_people_classify) <- name_list















